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ABSTRACT 


New  cloture  theorems  for  shock  models  In  reliability  theory  are  presented. 
If  the  niaber  of  shocks  to  failure  and  the  times  between  the  arrivals  of 
shocks  have  probability  distributions  of  phase  type,  then  so  has  the  time  to 
failure.  PH-distributions  are  highly  versatile  and  may  be  used  to  model 
many  qualitative  features  of  practical  interest.  They  are  also  well-suited 
for  algorithmic  implementation.  The  computational  aspects  of  our  results  are 
discussed  in  some  detail. 
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1.  INTRODUCTION 


Shock  models  which  relate  the  life  distribution  H(*)  of  a  device, 
subject  to  failure  by  shocks  occur lag  rand only  in  tine,  have  received  con¬ 
siderable  attention  In  recent  years.  If  Is  the  probability  that  the 

device  survives  k  >  0,  shocks  and  N(t)  Is  the  random  number  of  shocks 
la  (0,t],  the  survival  probability  H(»)  ■  1  -  H(*),  of  such  a  device  is 
given  by 

(1)  H(t)  -  -  Z  P.  P{N(t)-k)  . 

N(t)  k-0  k 

The  most  general  shock  models  are  those  that  correspond  to  (1) ,  such  that 
(N(t):  t  ^  0}  is  a  general  counting  process  and  1  >_  FQ  >_  >_  i>2  >_  . . .  . 

Interest  in  and  published -reaulM  for  shock  models  center. around  proving  , 
that,  subject  to  suitable  a as  motions  on  the  point  process  N(t)  of  shocks, 
various  reliability  characteristics  of  the  shock  resistance  probabilities 
Pk  are  inherited  by  the  survival  probability  H(»)  in  continuous  time. 

The  first  systematic  treatment  of  such  shock  models  was  given  by  Esary, 
Marshall  and  Prose han  [5],  when  N(t)  is  a  homogeneous  Poisson  process. 
A-Hameed  and  Proschan  considered  the  cases  when  N(t)  is  a  non -homogeneous 
Poisson  Process  [1]  and  a  non-statlonary  pure  birth  process  [2].  Block  and 
Savits  [4]  treated  the  case  when  the  lnterarrlval  time  between  shocks  is 
NBUE  (NBUE)  or  NBU(NWU)  and  Thall  (6]  derived  Interesting,  but  comparatively 
weaker,  results  idten  N(t)  is  a  clustered  Poisson  process. 

In  this  paper,  we  obtain  preservation  theorems  for  the  shock  model  (1) 
when  p^  is  of  phase-type  and  so  is  the  distribution  of  the  lnterarrlval 
time  between  shocks.  N(t)  is  then  a  phase  type  renewal  process  [7],  The 
relevance  of  phase  type  distributions  (henceforth  abbreviated  as  PH -distri¬ 
butions)  to  the  algorithmic  analysis  of  the  time  dependent  behavior  of 


of  stochastic  models  has  baan  discussed  by  Neuts  in  a  aeries  of  papers 
starting  with  [6].  A  comprehensive  treatment  may  be  found  in  Chapter  2  of 
[8].  PH-distributlons  provide  an  alternative  point  of  departure  in  modelling 
real  life  distributions  without  the  classic  memoryless  property  and  with 
possible  proper  unlmodallty  or  multimodality.  PH-distrlbutions  include  the 
exponential.  Erlang  and  hyper-exponential  distributions  as  very  special 
cases.  In  addition,  they  have  the  desirable  property  of  being  closed  under 
both  finite  convolutions  and  mixtures,  a  feature  possessed  by  none  of  the 
well-known  non-par ame trie  classes  of  life  distributions. 

In  Section  2,  the  basic  properties  of  PH-dlstributions,  needed  in  the 
sequel,  are  briefly  reviewed.  The  main  theoretical  results  are  discussed 

•  ••  •  •  «  s  *  •  •  •  • 

in  Section  3.  Algorithmic  considerations  are  presented  in  Section  4. 

2.  PH-DISTRIBUTIONS 

A  density  {pfc}  on  the  nonnegative  integers  is  of  phase  type  if  and  only 
if  there  exists  a  finite  Markov  chain  with  transition  probability  matrix  P 
of  order  r  +  1  of  the  form 

S  S' 

P  -  ~ 

0  1 

and  initial  probability  vector  [£,8^^] ,  such  that  {p^}  is  the  density 
of  the  time  till  absorption  in  the  state  r  +  1.  The  matrix  I  -  S  is 
nonsingular  and  the  stochastic  matrix  S  +  (1-8^)  _1S*  *j3  may  be  chosen  to 
be  irreducible. 

The  density  Cpfc>  is  given  by  pQ  -  8^,  and  Pk  ■  1  Sk_1  S* ,  for 
k  >_ 1.  In  this  paper  {p^}  will  be  the  density  of  the  number  of  shocks  to 
failure  in  a  reliability  shock  model.  Ve  will  assume  throughout  that 


3 


8 


nrt-1 


0 


We  also  clearly  have  chat 


—  V 

P.  ■  Z  p.  •  £  S  £,  for  k  >  0  . 

*  v-k+1  * 

The  aean  y^  of  {p^}  la  given  by  £(I-S)”*e  . 

A  probability  distribution  F(*)  on  [0,«)  is  of  phase  type  if  and  only 
if  there  exists  a  finite  Markov  process  with  generator  Q  of  the  form 


Q  - 


T 

0 


I* 

0 


with  initial  probability  vector  ] ,  such  that  F( •)  is  the  distri¬ 

bution  of  the  time  till  absorption  in  the  state  m  +  1.  The  matrix  T  is 
nonsingular  and  the  generator  T  +  (1-u^^)""1  T**o  may  be  chosen  to  be  irre¬ 
ducible.  The  distribution  F( •)  is  given  by 

(2)  F(x)  “1-u  exp  (Tx)  e_  ,  for  x  >.  0  . 

We  shall  denote  1  -  F(x)  by  F(x).  The  mean  A^  of  F(*)  is  given  by 

~  2.  T  1  «.•  Th®  pairs  <£,T)  and  (£,S)  are  called  the  representations 

of  F(«)  and  {p^}  respectively.  Renewal  processes  in  which  the  underlying 
distribution  F(*)  is  of  phase  type  were  discussed  in  [7]. 

Many  derivatives  related  to  PH-dlstrlbutlons  involve  the  Kronecker 
product  L  9  M  of  two  matrices  L  and  M.  This  is  the  matrix  made  up  of 
the  blocks  {L^M}.  Provided  the  matrix  products  are  defined,  we  have  that 

(3)  (L9M)  (K0R)  •  LK  0  MH  . 


This  property  is  repeatedly  used  in  the  sequel. 
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3.  CLOSURE  THEOREMS 

We  first  consider  the  Esary-Marshall-Proschan  (E.M.P.)  shock  model 
[3,5]  in  which  (N(t)}  is  a  Poisson  counting  process  of  rate  X. 

Theorem  1 

If  the  number  of  shocks  to  failure  has  a  discrete  PH-density  (p^.k  >_  0} 
with  representation  (S,S),  then  the  time  to  failure  in  the  E.M.P.  model  has 
a  continuous  PH-diatribution  H( •)  with  representation  [jB,X  (S-I)  ]  . 

Proof 

—  Jq 

Since  ■  6  S  e ,  for  k  0,  we  obtain 

QO 

S(t)  -  1  e“Xt  1  Sk  e  -  B,  exp  [X(S-I)t]  e  ,  for  t  >0  . 

k-0  s 

This  proves  the  stated  result. 

A  number  of  Interesting  quantities  may  now  be  expressed  in  computationally 
convenient  forms.  The  j-th  noncentral  moment  of  H(‘)  is  given  by 

(4)  uj  -  Jl  X"j  j3(I-S)-J  e  ,  for  ]  >  1  . 

The  density  h(t)  -  H' (t),  is  given  by 


(5)  h(t)  ■  X  ^  exp  [X CS— I) t ]  S_®  ,  for  t  >_  0  , 

_-l 

and  the  failure  rate  r(t)  *  h(t)  H  (t)  ,  equals 


(6) 


Bexp(XtS)S* 

"  x  5S5T5ts5T  •  for  'i* 


Theorem  1  is  a  particular  case  of  a  more  general  result  in  which  the 
arrivals  of  shocks  occur  according  to  a  PH-renewal  process  [7] .  This  result 
is  proved  next. 

Let  the  lnterarrlval  time  distribution  ?(•)  be  of  phase  type  with 


irreducible  representation  (a,,T)  of  order 


■  1  -  a  e ,  is 


positive,  a  geometrically  distributed  number  of  shocks  occur  slaultsneously 


5 


at  each  shock  epoch.  As  in  [7],  we  introduce  the  matrices  P(k,t),  k  >_0, 
t  >_  0,  which  satisfy  the  system  of  differential  equations 


(7) 


P*(0,t)  -  P(0,t)T  , 


v-1 


P'(k,t)  -  P(k,t)T  +  t  P(k-v,t)T#a  , 

v«l 


k  >  1  , 


for  t  >  0,  with  initial  conditions  P(k,0)  ■  I,  for  k  ^  0  . 

The  element  P^fk.t)  is  the  conditional  probability  that  the  Markov  pro¬ 
cess  with  generator  Q*  »  T  +  (1-a^^)  *  T*<x,  is  in  the  state  j  at  time 
t  and  that  k  shocks  have  occurred  in  (0,t],  given  that  it  started  in 


the  state  i  at  time  0. 

The  Markov  process  Q*  may  be  started  according  to  any  initial  pro¬ 
bability  vector  "i/  With  -  £—  (l-a^^T1.  a  ,  .the  ?H-reneval  process  is 
started  immediately  after  a  renewal  epoch.  With  **  -X^  ^oT  1,  where 
X^  *  -a  T_1  e,  is  the  mean  time  between  shocks,  we  obtain  the  stationary 
version  of  the  PH-renewal  process. 

Theorem  2 

If  the  shocks  occur  according  to  a  PH-renewal  process  with  underlying 
representation  <o,T)  and  the  process  Q*  is  started  according  to  the  pro¬ 
bability  vector  x.  an<*  1*  the  probability  density  {p„, )  is  of  phase  type 
with  representation  (8,S)  of  order  r,  then  the  distribution  H(»)  is  of 
phase  type  with  the  representation 

C8)  je,  -  X  •  A  » 

k  -  t  •  i  +  r®  •  S  • 


of  order  rm  . 


Proof 


By  Che  lav  of  CoCal  probability,  we  have 

00 

(9)  H(t)  -  X  1  P(k,t)e  •  B  Sk  e 

k-0 

m 

»  (X»6)  Z  P(k, t)  9  Sk  (eOe) 
k-0 

-  (x*D  z(t)  (*£«.)»  for  c  >  o  . 


The  matrix  z(t)  satisfies 


Z' (t)  -  E  P' (k,t)  9  Sk  -  E  P(k,t)  T  9  Sk 


k-0 


k-0 


+  E  E  o^1  P(k-v,t)  T*«  9  Sk 
k-1  v-1  ®fl 


Z(t)  (T9I) 


+  E 
k-0 


P(k,t)  T°a  9  Sk+1  (I-a^S)-1 


-  Z(t)  [T9 1  +  T°a  9  (I-a^jS)”1*]  , 
and  clearly  Z(0)  -  I  9  I  . 

This  implies  that  Z(t)  *  exp(Kt),  for  t  >  0.  Upon  substitution  into 

(9) ,  the  proof  is  complete. 

Particular  Cases 

1.  If  the  number  of  shocks  to  failure  is  geometrically  distributed,  i.e. 

—  1c 

Pfc  •  0  ,  for  k  >  0,  then 

k  _i 

(10)  H(t)  -x  z  P(k,t)flKe  -  x  < {T+(l-ea  ..)  •  fajt}  e  , 

k-0 

for  t  >0  . 

2.  In  the  maximum  shock  model,  failure  occurs  if  and  only  if  a  shock  occurs 
vhose  magnitude  exceeds  a  critical  randomized  threshold  Y  vlth  distribution 


G(*).  If  che  magnitudes  of  successive  shocks  are  independent  with  common 
distribution  F(»),  then 


(11)  Pk  -  j  F^x)  dG(x)  ,  for  k  >  0  . 

0 

It  follows  from  (10)  that 

(12)  H(t)  -  J  X  exP  ttT  +  (l-am+>1F(x))“1F(x)T0a]t}  e  dG(x)  , 

0 

for  t  >_Q,  so  that  H(«)  is  a  mixture  of  PH-distributions .  If  G(')  is 
a  discrete  distribution  with  finite  support,  then  H(*)  itself  is  of  phase 
type. 

3.  In  the  cumulative  damage  model,  the  damages  are  additive.  With  the  same 

distributions  F(«)  and  G(«)  as  in  the  preceding  model,  we  obtain 

00 

(13)  Pk  -  |  F<k) (x)  dG(x)  ,  for  k  >  0  . 

0 

If  the  distribution  G(«)  is  of  phase  type  with  representation  d,L),  then 

00 

Pk  •  |  G(x)  dF(k)  (x)  -  E  G(xx  +  ...  +  x^ 

0 

■  E  !  exp  [L(x^  +  . . .  +  Xj^)  ]e_  “  Ak  e_  , 

00 

where  A  ■  j  exp  (Lx)  dF(x) .  It  is  readily  seen  that  A  is  a  substochastic 
matrix  of  spectral  radius  less  than  one.  The  density  {pk>  is  therefore  of 
phase  type.  If  the  shocks  occur  according  to  a  PH-renewal  process.  Theorem 
2  may  be  applied  to  evaluate  H(t).  The  matrix  A  is  obtained  by  numerical 
Integration  for  general  distributions  F(')>  If  F(»)  itself  is  of  phase 
type  with  representation  (o_, R) ,  then 
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(14)  A  ■  j  exp  (Lx)o  exp(Rx)R*dx  ■  (I8o)  J  exp(Lx)9exp(Rx)  dx  ( 10  R° ) 

0  0 
,  -  -d ®o)  [Lei  +  ieR]-1  (ieR#)  . 

The  eigenvalues  of  L  and  R  all  lie  in  Che  open  left  half-plane.  The 
same  Chen  holds  crue  for  Che  Kronecker  sum  Lei  +  I0R,  so  chac  Che  Inverse 
exiscs. 

The  nonnegaclve  recCangular  macrix  V  ■  -(Lei  +  I0R)-1  (iejR0),  may 
easily  be  compuced  by  solving  Che  sysCem 

(Lei  +  ieR)v  -  -  ieR° 
by  block  Gauss-Seldel  iteration. 

4.  ALGORITHMIC  ASPECTS 

We  shall  discuss  Che  compuCadon  of  Che  funcdon  H(t),  which  is  given 
by  Theorem  2.  Ic  readily  follows  from  (1)  chac  che  mean  h^  of  H(*)  is 
given  by  XJp^,  where  A|  and  are  Che  means  of  {p^}  and  F(*)  re¬ 

spectively,  whenever  Che  PH-renewal  process  of  arrivals  is  scarred  aC  a 
renewal  epoch.  Wlch  general  lniclal  condicions,  Che  mean  hj^  is  given  by 
X’y^  +  X^  -  X^,  where  X^  -  -  ^  T~*  e  . 

Knowledge  of  che  mean  h^  of  H(«)  is  useful  in  decermining  che 
incerval  over  which  we  wish  Co  evaluaCe  H(c) .  We  may  e.g.  wish  co  choose 
che  mean  as  a  convenienC  unit  of  time.  This  is  accomplished  by  replacing 
K  by  hjK.  A  different  rescaling  may  be  chosen  if  Che  elements  of  hjK 
are  very  large  or  if  a  different  time  scale  is  desirable  for  che  practical 
problem  at  hand. 

We  now  assume  chac  Che  matrix  K  has  been  appropriately  rescaled.  The 
function  H(t)  is  compuced  by  numerical  integration  of  che  system  of  linear 


differential  equations 


(15) 


v'  (t)  •  v(t)  K  >  for  t  >_  0  , 

v(0)  -  1  8  £  , 


and  setting  H(t)  ■  v(t)£,  for  t  £  0  . 

It  is  convenient  to  partition  the  vector  v(t)  as  [v.(t),  ....  v  (t)], 

—  “"l  m 

where  the  vectors  v^(t)  are  r-vectors.  We  also  set  M  •  (I-am| ^S)  ^S. 


The  system  (IS)  may  then  be  rewritten  as 


m 

,  E  v  (t)  T°  M  , 
1  .  —v  v 

biri  _ 


for  1  <_  j  <_  m.  This  system  may  be  conveniently  solved  by  a  classical  inte- 
■  •  •  •  •  •••  4 

gration  procedure,  such  as  Runge-Kutta.  We  see  that  the  vector 

I  m 

I  v  (t)  T°  M  does  not  depend  on  j  and  needs  to  be  evaluated  only  once 

hr1  11 

in  each  computation  of  the  right  hand  sides  of  (16) . 

In  many  PH-distributlons  of  practical  interest,  such  as  e.g.  finite  mix¬ 
tures  of  Erlang  distributions,  the  order  m  of  T  may  be  large,  but  T, 

T°  and  have  very  few  nonzero  entries.  It  is  then  advantageous  to  write 
a  special  purpose  subroutine  to  evaluate  the  right  hand  side  of  (16).  By 

so  exploiting  the  sparsity  of  T,  T#  and  o,  it  is  possible  to  reduce  the 

computation  time  greatly.  The  mean  h’,  or  in  general  the  scaling  factor 
used  in  selecting  the  time  unit,  may  also  be  utilized  to  choose  the  step  size 
h  in  the  numerical  integration  of  the  system  (16).  In  similar  problems,  we 
have  usually  made  two  runs  at  least,  one  with  1/50  of  the  time  unit  and 

one  with  1/100  of  the  time  unit.  If  the  results  at  corresponding  time 

points  are  not  sufficiently  close,  further  runs  with  smaller  steps  are  made. 
The  computation  times  of  such  runs  increase  rapidly  and  efficient  programming 


is  desirable.  Ocher  methods  with  a  variable  step  size  and  error  control 
may  also  be  Implemented.  These  classical  topics  in  the  numerical  integration 
of  ordinary  differential  equations  need  not  be  belabored  here.  In  all  cases, 
the  use  of  the  particular  structure  of  the  matrix  K  is  fully  worthy  of  the 
additional  programing  effort. 
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